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SUMMARY 

A time-marching finite-difference method has been used to solve the 
compressible Navier-Stokes equations for the three-dimensional wing-leading- 
edge shock impingement problem. The bow shock was treated as a discontinuity 
across which the exact shock jump conditions were applied. All interior 
shock layer detail such as shear layers, shock waves, jets, and the wall 
boundary layer are automatically captured in the solution. The impinging 
shock was introduced by discontinuously changing the freestream conditions 
across the intersection line at the bow shock. A special storage- saving 
procedure for sweeping through the finite-difference mesh has been developed 
which reduces the required amount of computer storage by at least a factor 
of two without sacrificing the execution time. Numerical results are 
presented for infinite cylinder blunt body cases as well as the three- 
dimensional shock impingement case. The numerical results are compared 
with existing experimental and theoretical results. 
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I. INTRODUCTION 

When an extraneous shock strikes the bow shock of a blunt body, a 
complicated flow pattern results which may contain imbedded shock waves, 
shear layers, expansion fans or supersonic jets. In addition, large in- 
creases in wall heat transfer and pressure can occur and have been mea- 
sured experimentally by many investigators. Impinging shock flow fields 
may occur on hypersonic aircraft or reentry vehicles such as the Space 
Shuttle. 


A. Classification of Shock Impingement 

The different flow fields which can result have been categorized by 
Edney [l] into six basic types of shock interference patterns. These 
patterns are shown in Figures 1-6. Impinging shock strength, freestream 
Mach number, body shape, and position of impingement (with respect to the 
subsonic flow region on the blunt body) are the primary factors in deter- 
mining the type of shock impingement. 

A Type I interference pattern (.shown in Figure 1) occurs when the 
impinging shock strikes the bow shock well above the upper sonic line. 

In this case the impinging and bow shocks are of opposite families. The 
shear layer emanating from the intersection point does not strike the 

i 

body and therefore, does not contribute to increases in wall pressure or 
local heating. However, the transmitted impinging shock can strike the 
body. If strong enough, this transmitted shock can cause boundary layer 
separation or transition from laminar to turbulent flow. Increases in 
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Figure 1. Type I shock interference 
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Figure 2. Type II shock interference 
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Figure 4. Type TV shock interference 




7 



BOW SHOCK 


Figure 6. Type VI shock interference- 
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wall pressure and heat transfer can then result from this shock-houndary 
layer interaction. Heating amplifications as large as ten can result 
from a Type I interference. 

A Type II interference pattern (shown in Figure 2) occurs when the 
impinging shock strikes the bow shock just above the upper sonic line. 

In this case the impinging and bow shocks are also of opposite families. 
The transmitted impinging shock strikes the body surface resulting in 
possible boundary layer separation or transition. In addition, two shear 
layers are formed which may strike the body to cause additional heating 
and pressure peaks. The shear layer emanating from point B is very weak 
and can cause only small heating and pressure amplifications. On the 
other hand, the shear layer emanating from point A is strong and can 
cause large increases in heating and pressure. Heating amplifications as 
large as five can be expected from a Type II interference. 

A Type III interference pattern (shown in Figure 3) occurs when the 
impinging shock strikes the bow shock just inside the upper sonic line. 

A shear layer emanates from the intersection point and strikes the body 
causing large increases in heat transfer and wall pressure at this point. 
This shear layer separates regions of subsonic and supersonic flow. The 
supersonic flow is turned parallel to the body by an oblique shock which 
exists between the bow shock and the body. Heating amplifications as 
large as ten can be expected for a Type III interference. 

When flow conditions for the Type III case become such that the 
oblique shock no longer is able to turn the flow parallel to the body, 
a Type TV shock interference is formed and is shown in Figure 4. This 



9 


interference pattern occurs when the impinging shock strikes the near- 
normal part of the bow shock. Imbedded in the subsonic portion of the 
blunt body flow is a supersonic jet which is terminated by a normal jet 
bow shock just above the body. The Type IV interference pattern pro- 
duces the largest heating and pressure peaks due to the jet impingement 
upon the body. Heating peaks as high as 17 times the undisturbed stag- 
nation value and pressure peaks as high as 8 times the undisturbed stag- 
nation value have been measured [2], 

A Type V interference pattern (shown in Figure 5) occurs when the 
impinging shock strikes the bow shock just below the lower sonic line. 

This case is nearly a mirror image of the Type II case, except the imping- 
ing and bow shocks are of the same family and a thin supersonic jet eman- 
ates from the impingement point (point B) instead of a shear layer. Heat- 
ing amplifications as large as five can be expected for a Type V inter- 
ference. 

A Type VI interference pattern (shown in Figure 6) occurs when the 
impinging shock strikes the bow shock well below the lower sonic line. 

This interaction is analogous to the Type I case except that the impinging 
and bow shocks are of the same family and an expansion fan is transmitted 
instead of a shock. The interaction of the expansion fan and the wall 
boundary layer results in a small decrease in the wall pressure and heat 
transfer. 

Further details about the different shock interference patterns are 
given in the original work by Edney. In addition, an informative discus- 
sion on shock interference classification is given by Keyes and Hains [2], 
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A two- dimensional shock interference geometry is shown in Figure 7a. 
This geometry can produce any of the six interference patterns. The im- 
pinging shock is planar and oriented such that the intersection line is 
parallel to the axis of the cylinder, causing the flow in each z-plane to 
be identical. This geometry is important because it qualitatively models 
the flow field in many three-dimensional cases. 

A more physically realistic geometry occurs when a planar impinging 
shock intersects a circular cylinder as shown in Figure 7b. This con- 
figuration can occur when the bow shock from the nose of a vehicle strikes 
the wing leading edge bow shock. A fully three-dimensional calculation 
must be performed to obtain this flow field. The sweep angle of the cyl- 
inder is a critical parameter in this case because it controls the inter- 
ference pattern type. For zero or small sweep angles a Type IV interfer- 
ence pattern will result. For moderate sweep angles (around 30°) a Type V 
interference pattern will occur. For larger sweep angles a Type VI inter- 
ference pattern will result. Because swept forward wings generally do 
not exist. Types I, II, and III are not considered for the wing leading 
edge case. Further discussion about this geometry and some of its flow 
field detail is presented in the following sections. 

B. Literature Review 
1. Experimental shock impingement studies 

The large number of publications appearing in the literature on 
shock interference is indicative of the amount of work being done in this 
area. Because of the complexity of the problem, most studies have been 
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experimental or have used a highly simplified model in a theoretical 
analysis. 

Edney [l], Keyes and Hains [2], and Korkegi [3] provide an extensive 
survey of the literature on the shock interference problem. References 
Cl] and [4-29], which are reviewed by Keyes and Hains, represent the bulk 
of the experimental effort from 1961 to 1972. Many different experi- 
mental approaches have been used to study the problem of shock interfer- 
ence. The hemisphere and leading edge impingement geometries have re- 
ceived most of the attention. The effects of Reynolds number, Mach num- 
ber, impinging shock strength, and sweep angle have been recorded through 
various pressure, heat transfer, and flow visualization measuring tech- 
niques. The scatter in the experimental data is unusually large, espe- 
cially the peak values of heat transfer and pressure. This can be attrib- 
uted to the very small region over which an impingement interaction 
occurs, making accurate measurements difficult. 

As described earlier, Edney in 1968 published a classical paper on 
shock impingement in which six basic shock interference patterns were 
classified (see Figures 1 - 6 ). In addition, new techniques for heat trans ‘ 
fer measurement were developed which improved the accuracy of the peak 
heating measurements. In addition to heat transfer measurements, pressure 
measurements and Schlieren photographs were obtained for hemispheres, 
leading edge fins, wedges, and flat-faced cylinders. 

Keyes «nd Hains [2] and Hains and Keyes [30] published the results 
of extensive experimental investigations of shock impingement. All six 
impingement types due to Edney were investigated. Planar shocks were 
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allowed to impinge upon the bow shocks surrounding both hemispherical 
and swept fin models. Heat transfer measurements were made by using a 
phase-change coating technique. This technique has the {•bii.ity to detect 
heating peaks over very narrow regions. The heat transfe* results, 
however, are subject to numerous and often large sources of error associ- 
ated with the measurement of 1) thermophysical properties of the model 
material, 2) melt temperature, 3) melt time, 4) initial model temper- 
ature, and 5) adiabatic wall temperature distribution. Peak heat trans- 
fers of 17 times the undisturbed stagnation values and peak pressures of 
8 times the undisturbed stagnation pressures were measured for the most 
critical (Type IV) hemispherical case. In an additional report by Keyes 
[31], off -center-line heating rates were measured for both hemispheres 
and swept fins. Heating* rates higher than the undisturbed stagnation 
value due to shock impingement were found over significant portions of 
the off-center-line body surface. 

Doughty et al. [32] determined angle of attack effects on shock 
impingement heating rates by testing several delta-wing space shuttle 
models. The same phase-change 'coating technique was used in this inves- 
tigation as was used in references [2] and [30] , At low angles of attack 
(20°) , the heating rate increases (due to shock impingement) were re- 
stricted to the vicinity of the wing leading edge. Heating rates on the 
windward side of the trailing edge were almost unaffected by the shock 
impingement. At moderate angles of attack (40°) increased heating rates 
due to shock impingement were experienced over a very large portion of 
the windward side of the wing. No mention was made of the heating rates 
on the leeward side of the wing. 
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Other recent experimental studies were performed by Bertin and 
Hinkle [33], Palko [34], and Craig and Ortwerth [35], Bertin and Hinkle 
used a double wedge model to approximate the wing leading edge problem. 
Palko conducted an extensive series of experiments also on the wing lead- 
ing edge problem. Stagnation line pressure and heat transfer measurements 
and Schlieren photographs were obtained for systematic variations of 
freestream Mach number, leading edge sweep angle, and impinging ';hock 
strength. 

Craig and Ortwerth investigated the leading edge shock impingement 
problem hut with a different impinging shock orientation. The orientation 
studied was that of the two-dimensional type shown in Figure 7a. This 
orientation exactly duplicates the type of flow fields shown in Figures 
1-6. Pressure and heat transfer measurements and Schlieren photographs 
were obtained. 

2. Theoretical shock impingement studies 

Many authors from references [4-29] have offered theoretical models 
and methods of peak heating prediction. Most of these authors, however, 
do not present a conclusive understanding of the shock impingement prob- 
lem. The results, therefore, are in most cases not satisfactory [2J. 

Edney established flow models for each of the six shock interference 
patterns. By determining the nature of the shock interference flow- 
fields, Edney established the mechanisms which produce the pressure and 
heating peaks. These peaks are caused by one or more of the following 
phenomena: 1) shock -boundary layer interaction, 2) free shear layer 

attachment, and 3) supersonic jet impingement. From these flow models 
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semi empirical methods for predicting the approximate flow fields were 
developed. Edney presented two methods, a graphical method and a numeri- 
cal method. 

The graphical approach, referred to as the heart diagram method, is 
based on a plot of pressure versus flow deflection angle and can be used 
to describe the flow fields for each of the six shock interaction types. 

The numerical method presented by Edney utilizes the oblique shock 
relations and an iterative procedure. Types I and VI can be solved 
exactly while Types II through V require semiempirical values for certain 
characteristic lengths. 

Crawford [36] improved Edney* s graphical approach by removing the 
iteration procedure. This was accomplished by plotting a family of 
pressure deflection curves with the pressure ratio on a logarithmic 
scale and flow deflection on a linear scale. 

Bramlette [37] improved Edney* s numerical scheme with two simplify- 
ing assumptions. This modified technique can only be applied to Type III 
and Type IV interactions, but removes the iterative portion of the solu- 
tion without sacrificing the accuracy. 

Many authors have presented correlations for peak heating due to 
shock-boundary layer interactions and free shear layer attachment [38-433. 
The most useful shock -boundary layer correlation, developed by Markarian 
[38] and used for shock impingement problems by Keyes and Morris [43], is 
based upon the inviscid pressure rise across the interaction region. The 
basic expression is given by 
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where and are the peak and undisturbed values of the heat trans- 
fer coefficient, respectively, p p and are the peak and undisturbed 
values of wall pressure, respectively, and C and m are empirical con- 
stants which depend upon the nature of the boundary layer (i.e., lami- 
nar, transitional, or turbulent) . 

The boundary layer reattachment correlation for peak heating, pre- 
sented by Bushnell and Weinstein [42] , was used by Keyes and Morris for 
the shock impingement free shear layer attachment problem, 

Keyes and Hains [ 2 ] and Hains and Keyes [3fl] present the results of 
an extensive theoretical study of all six types of shock impingement. 

The flow models introduced by Edney and the peak heating correlations 
described above have been utilized and extensively compared Tfith experi- 
mental results. These theoretical methods for each impingement type have 
been coded into computer programs and are listed along with detailed 
descriptions in reference [44] . These theoretical methods provide fast, 
easily obtainable, and in most cases (where real gas effects are negligi- 
ble and good flow visualization photographs are available) adequate an- 
swers for the peak heating and pressure values resulting from the shock 
impingement. However, many disadvantages for these theoretical methods 
exist [2l: 1) The models used are local, two-dimensional models and do 

not include three-dimensional effects. 2) The calculations are based 
upon flow visualization photographs or other empirical sources. 3) 

Shear layer growth or oblique jet impingement are not included. 
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4) Gas chemistry is not included. 

Bertin et al. [45] presented a theoretical study of shock impinge- 
ment (Types V and VI) on a double wedge model. Real gas effects were 
included. The method utilized oblique shock relations. Included were 
perfect gas results, variable gamma (Y) results and results for air in 
chemical equilibrium in which all properties were obtained by a table 
look-up scheme. Because of real gas effects, Bertin et al. [45] con- 
cluded that extrapolation of wind tunnel data to flight conditions may 
yield inaccurate results (e.g., in the case of the Space Shuttle). 

3. Numerical shock impingement studies 

The numerical schemes for computing blunt body flow fields have had 
a pronounced effect upon shock impingement numerical solutions. There- 
fore, a brief discussion of blunt body numerical solutions is presented 
first. 

Many schemes have been developed to solve the classical blunt body 
problem in a supersonic flow. Discussion and classical theory on the 
blunt body problem can be found in Hayes and Probstein [46J . Daywitf and 
Anderson [47]present a review of blunt body solutions in which special 
attention is given to the time-dependent, finite-difference approach. 

Two basic variations utilizing the time -dependent approach have been 
developed. The simplest is the so-called ’'shock-capturing” technique 
[48-50] . In this technique the detached bow shock is automatically com- 
puted with no a priori knowledge of Its location or even existence. 

The later technique is the so-called ’’shock-fitting" technique [5l] . 
The bow shock in this approach is assumed to ba a boundary of the 
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computational mesh across which the exact shock jump conditions 
(Rankine-Hugoniot relations) are applied. 

Another study [523 presented a viscous blunt body solution for a 
gas in chemical nonequilibrium. The time-dependent shock fitting 
approach was used in conjunction with the Navier-Stokes equations. In 
this approach, not only the inviscid flow field, but also the boundary 
layer and the chemical species concentrations were computed. 

Other numerical techniques for blunt body flows exist, but are not 
discussed here. Instead, the topic is changed to blunt body flows with 
an impinging shock. 

Tannehill and Holst [53] and Tannehill et al. [54] both obtained 
two-dimensional shock impingement solutions using the time -dependent 
approach. In both papers, the complete Navier-Stokes equations were 
solved. The geometry used is shown in Figure 7a. 

The first paper obtained a low Reynolds number solution simulating 
high altitude conditions by using the shock capturing approach. In addi- 
tion, a preliminary study of fitting the bow shock and capturing all 
shock layer detail was presented. 

The second paper presented several moderate Reynolds number solutions 
utilizing the shock fitting method. The impinging shock was introduced 
by discontinuous ly changing the freestream conditions at the bow shock. 

Holst et al. [SS] presents a brief discussion of the present results. 

This completes the literature review. A brief presentation of the 
methodology used in the present study for computing shock impingement 
flow fields is presented in the next section. 
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C. Present Methodology 

In the present study the three-dimensional wing leading edge shock 
impingement flow field is numerically computed by using a time -dependent 
finite-difference method. The choice of a numerical method to solve shock 
impingement flow fields seems well suited because of the complicated type 
of problem involved. 

The three-dimensional formulation is important. The inclusion of the 
third dimension provides a "relief effect" which can drastically alter 
too -dimensional results. 

The time-dependent approach was chosen because a subsonic region may 
exist along the stagnation plane of the wing leading edge. Even if the 
cross-flox? velocity is supersonic before impingement, a subsonic region 
may form after the impinging shock has been introduced. Since the gov- 
erning time-dependent equations remain a hyperbolic -parabolic set for 
both subsonic and supersonic flows, all cases can be solved as an initial - 
value problem where the steady-state solution is approached asymptotically 
with time. 

A standard approach for the numerical solution of a blunt body prob- 
lem is to first compute an inviscid solution using the Euler equations 
and then with the resulting pressure distribution compute a boundary layer 
solution. This approach is invalid for the present case because of the 
various types of boundary layer interactions caused by the shock impinge- 
ment. The Navier-Stokes equations, however, are valid in such situations 
and have been chosen in the present study. In particular, the laminar, 
compressible form of the Navier-Stokes equations have been chosen. 
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Because of the possible existence of shock waves in the computational 
domain, a conservative form of the governing equations has been used. 

The bow shock has been treated as a boundary of the computational 
domain across which the exact shock jump conditions (Rankine-Hugoniot 
relations) have been applied. This is the so-called "shock fitting" 
procedure. All interior detail such as shear layers, shock waves, super- 
sonic jets, and the wall boundary layer are automatically captured. This 
is the so-called "shock capturing" procedure. The bow shock is strong 

4 

enough to make an entirely captured solution impractical. On the other 
hand, fitting any additional Bhock layer detail (such as the transmitted 
shock) is generally too complicated (because of the transmitted shock- 
boundary layer interaction). Therefore, fitting the bow shock and cap- 
turing all shock layer detail seems to be a good approach for solving 
the problem at hand. 

For the shock impingement problem, a distinct advantage for numerical 
methods exists over theoretical or experimental methods. Scale effects 
must be considered when extrapolating experimental data to flight condi- 
tions. Present theoretical methods require empirical inputs. However, 
numerical methods are able to predict flow fields for flight conditions 
without empirical inputs. This is a particular advantage for numerical 
methods when real gas effects enter the problem. 

Descriptions of the governing equations, the numerical method, and 
the numerical results are given in the following sections of this report. 
The results consist of a series of blunt body infinite cylinder solutions 
and a three-dimensional wing leading edge shock impingement solution. 
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The present results are compared with theory and experimental results 
when possible. 
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II, GOVERNING EQUATIONS 


A. Navier-Stokes Equations in Cylindrical Coordinates 


The equations governing the flow of a compressible, viscous fluid 
in the absence of body forces and electromagnetic effects, and written 
in weak conservation- law form using three-dimensional cylindrical coor- 
dinates. are given by [56,57] 
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The Navier-Stokes expressions for the components of the stress 
tensor (*L .) and heat flux vector (q^) have been used in this study and 
are given by [56,58] 
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To complete this set of equations two perfect gas equations of 
state are used. In addition, Sutherland's equation and a constant 
Prandtl number assumption are used to compute coefficients of viscosity 
(p.)and thermal conductivity (k) , The resulting equations are given by 


p = 

pe( y- 1) 

2,11 

p “ 

pgRT 

2.12 


3/2 

C 1 T /<C 2 + T) 

2.13 

k = 

0 ^ [VPr 

2.14 


where for air in the perfect gas range the constants have the following 
values: 

gR = 287.024 m 2 /sec 2 - °K 
c - 1.2585 x 10 8 N-sec/m 2 - °K l/2 

c 2 =110.4°K 2.15 

c = 1004.585 m 2 /sec 2 - °k 
P 

Pr = 0.72 


B. Independent Variable Transformation 

Equations 2.1 through 2.10 are transformed from the physical domain 
(r, 0, z, t) into the computational domain ( rt, £, t) using the following 
independent variable transformation: 
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The derivatives associated with this transformation are given by 
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where r is the bow shock radius, r, is the body radius, and r q 
s b , 9 z 

are the shock slopes in the R- and z-directions, respectively, r S( _ is 
the radial shock velocity, and x is defined by 


r - r 



This transformation maps the z-plane between the bow shock and the 
blunt body into a rectangular region, stretches the radial distribution 
of grid points according to the function f, and stretches the axial dis- 
tribution of grid points according to the function g. 



The function f chosen for all cases considered here is given by 
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where 

c = ln (frr) 2 - 20 

The first and second derivatives of the radial stretching function 
(f) are given by; 


af _ 2g/c 
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Equation 2.19 refines the grid near the body and thus permits 
better boundary layer resolution. Figure 8 is a plot of the f stretch- 
ing function for various values of 0 which controls the amount of re- 
finement. The practical range for 0 is between 1 and 2 with smaller 
values giving larger amounts of refinement. When 0 is set equal to 
zero no radial stretching of the mesh occurs. 

The function g chosen for all cases considered here is simply that 
of a linear distribution (no stretching of the mesh) and is given by 


2.22 
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Figure 8. f stretching function distribution. 
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= 1.0 2.24 

dz 

,2 

= 0.0 2.25 

Sz 

The physical domain of the problem is shown in Figure 9a without 
an impinging shock for clarity. The computational domain which results 
after the independent variable transformation has been applied is shown 
in Figure 9b. The grid distribution appearing in a typical 0-plane in 
the physical domain is shown in Figure 10a. The grid distribution in 
a typical r|-plane, which is the computational counterpart of the physi- 
cal 0-plane, is shown in Figure 10b. These two figures illustrate what 
the independent variable transformation does to the physical domain. 

No matter how the physical domain changes or alters its shape during 
the computation, the computational domain retains its rectangular box 
shape with uniform mesh spacing throughout. 


C. Governing Equations in Computational Domain 


The 

variable 


final form of the governing equations after the independent 
transformation has been applied is given by 
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where 
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IMPINGING SHOCK 



b. TYPICAL n -PLANE (COMPUTATIONAL DOMAIN) 


Figure 10, A planar comparison of the physical and computational domains. 
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In these equations the radial position (r) is given by 
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**1 

where the function f (5) represents the inverse of f. For the present 
f (Equation 2.19) the inverse is given by 
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In addition Equations 2.7 and 2,11 through 2.15 must be included 
and are unaltered by the transformation. 

Equations 2.26 through 2.34 represent the full three-dimensional 
Navier-Stokes equations valid for compressible, laminar, Lime -dependent 
flows. These equations are solved by a numerical finite-difference 
procedure which is presented in the next section. 
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III. NUMERICAL METHOD 


A. Finite-Difference Scheme- 


Equation 2.25 is solved by MacCormack' s explicit finite-difference 
method [60-], which is composed of the predictor-corrector sequence given 
in Equations 3.1 and 3.2. 


predictor step : 


n+T n 
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corrector step: 
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where F i,i,k = F ( U i,j,k)’ F i!j!k “ F ( U i J,k)’ etc ‘ 1116 1>J * k subscripti 

and the n superscript correspond to the discretized position in the 


finite-difference mesh in the f-, n-, and T-directions such that 
§=iA£, rt= jArv C~ k ACi and T- nAT. A bar above the n superscript in- 
dicates a predicted value. Forward differences are used in the predic- 
tor step to approximate the SF/^g, AG/St^ and SH/dC derivatives vfaile 
backward differences are used for these terms in the corrector step. 


In the predictor step backward differences are used to approximate the 


37 


shear stress and heat flux derivatives appearing in the F, G, H, and D 
components while in the corrector step forward differences are used* 

This differencing procedure maintains second-order accuracy in both 
space and time [61]. 

The computational module for MacCormack’s method in three spacial 
dimensions is shown in Figure 11. Both the predictor and corrector 
step modules are shown. The symbol o is used to indicate mesh points 
at which values for F, G, H, or D are required and are used in differ- 
ences approximating the SF/S^j dG/Sri, and derivatives. The sym- 

bol -f is used to indicate mesh points requiring values for u , u , u , 

3T 0 Z 

or T and are used in differences approximating the shear stress and 
heat flux derivatives. Note that in both the predictor step (Equation 
3.1 and Figure 11a) and the corrector step (Equation 3.2 and Figure 11b) 
information from only three k-planes is required (k-l,k,k+l). In the 
predictor step (Equation 3.1) the subscripts k and k+1 appear explic- 
itly while the subscript k - 1 is implicitly imbedded in the shear stress 
and heat flux derivatives. In the corrector step (Equation 3.2) the 
subscripts k - 1 and k appear explicitly while the subscript k+ 1 is im- 
plicitly imbedded in the shear stress and heat flux derivatives. 

For the present study MacCormack's method is applied to the finite- 


difference mesh in four steps: 1) partial-predictor step, 2) completor- 

predictor step, 3) partial-corrector step, and 4) completor-corrector 
step. These steps are shown in Equations 3. 3-3.6. 



b. CORRECTOR STEP 


Figure 11. Computational modi 
spacial directions 
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partial-predictor step : 
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partial-corrector step: 
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where the pp superscript designates the partial-predictor step result 
and the pc superscript designates the partial-corrector step result. 

All other notation is identical with the notation used in the standard 
MacCormack method. Note that when steps 1 and 2 (Equations 3.3 and 3.4) 

it 

are combined the standard MacCormack predictor step (Equation 3.1) is 
obt$?,v»ed and that when steps 3 and 4 (Equations 3.5 and 3.6) are com- 
bined the standard MacCormack corrector step (Equation 3.2) is obtained. 
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After each completed predictor or corrector step the U vector 

■“i 
U 2 
U 3 
°4 
. U 5 

is decoded to obtain the primitive variables (p, u r » u Q > u z » e ) in the 
following manner: 

u r * V°1 

= V°i 

“z * V U 1 3,8 

e - Oj/Uj - (>Y + u 9 + u ^ 2 

P" V r(r s' r b> 

In all four steps shown in Equations 3.3 -3.6 information from only 
two k-planes is required in any single step. This aspect of the modi- 
fied MacCormack method allows for a unique, storage- saving sweeping pro- 
cedure. Ibis procedure is described in detail in the next section. 

B. Application of the Finite-Difference Scheme 
(Sweeping Procedure) 




MacCormack* s method is usually coded to compute predicted values 
at each grid point over the entire finite-difference mesh before the 
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corrector step calculation is initiated. This requires storage for 
both the predicted and corrected solution arrays during the computa- 
tional process. In addition, many codes store other intermediate re- 
sult arrays over the entire finite-difference mesh. For two-dimen- 
sional problems this extra storage may not be a critical restriction. 

For three-dimensional problems these extra arrays may require large 
quantities of bulk storage causing significant penalties in program 
efficiency. 

The storage restriction associated with MacCormack' s method for 
three-dimensional problems has been to a large extent alleviated in 
the present study. A special technique for sweeping through the 
finite -difference mesh has been developed which requires storage for 
only a single three-dimensional solution array. Besides the single 
three-dimentional solution array several smaller too -dimensional arrays 
(k-planes) are required to store intermediate results. This special 
sweeping procedure is shown in Figure 12 where it is displayed in five 
separate diagrams (steps 1-5). For this presentation it is assumed 
that the computer memory is divided into two types: 1) a larger, 

slower access large core memory (LCM) , and 2) a smaller, faster access 
small core memory (SCM) . Data exchange between small and large core 
memory is efficiently achieved by using block transfer. This is the 
basic architecture of many computing systems such as the GDC 7600, CDC 
STAR and Burroughs ILLIAC IV. The three-dimensional solution array re- 
sides in LCM and appears at the top of each of the diagrams in Figure 12. 
The intermediate result arrays reside in SCM and appear below the three- 
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Figure 12. Diagram of the three-dimensional sweeping procedure. 
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Figure 12. (Continued). 
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Figure 12. (Continued). 


45 


dimensional solution array in each of the diagrams. These arrays are 

indicated as planes in Figure 12 with operations between them being 

indicated by arrows. Each plane represents a set of five arrays, one 

array for each of the primitive variables ( p, u , u 0 , u , e) . The basic 

i r y 2 

sweeping direction is the k-direction (z-direction) . 

The first step of the sweeping procedure is shown in Figure 12a. 
First, two planes of the n level solution (a total of ten arrays) are 
read into SCM from the k = 1 and 2 planes of LCM. From the data in 
these two planes the k = 2 plane of partial-predictor values is com- 
puted. A secondary sweeping procedure is used to obtain these partial- 
predictor values. This sweeping procedure starts at the bow shock (i=l) 
and sweeps through the lc-plane of mesh points toward the body (i=NI). 

The purpose of this planar sweeping procedure is to reduce the number 
of k-plane intermediate result arrays. Boundary conditions for the 
n-fl and the n+1 level solutions at k - 1 are also applied in this 
step (see Figure 12a) . 

The next step continues the start up procedure and is shown in 
Figure 12b. Since the n level plane at k = 1 is no longer needed, it 
is overwritten with the n level plane at k = 2. Next, a new n level 
plane from LCM at k = 3 is read into SCM, keeping the number of n level 
planes in SCM at two. A plane of partial-predictor values at k = 3 and 
a plane of completor-predictor values at k = 2 can now be computed. 

The completor-predictor step at k = 2 uses the partial-predictor values 
which were computed and saved from the previous plane at k = 2. Then, 
using the k = l and 2 planes of n+1 level values and the k-2 plane of 
n level values, a plane of partial-corrector values at k = 2 is computed. 
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Finally the n+1 level plane at k=l is stored in the three-dimensional 
solution array by overwriting the k = l plane of the n level solution. 
Planes containing the n+1 level solution in the three-dimensional solu' 
tion array in the diagrams of Figure 12 are shaded while n level planes 
remain unshaded. 

A recursive process for the general kth plane has the following 
steps (see Figure 12c): 1) read the k+2 plane of the n level solution 
from LCM into SCM, 2) compute a plane of partial-predictor values at 
k+2, 3) compute a plane of comp le tor -predictor values at k+1, 4) com- 
pute a plane of partial-corrector values at k+1, 5) compute a plane 
of completor-corrector values at k, and 6) store the kth plane of the 
n+1 level solution in the three-dimensional solution array by over- 
writing the kth plane of n level solution. The recursive process 
begins at k = 2 and continues until k = NK-3 where NK represents the 
maximum value of k. Note that by constantly overwriting unnecessary 
planes of storage the maximum number of intermediate result planes in 
SCM never exceeds seven; two for the n level solution, three for the 
n + 1 level solution and two for the n+1 level solution. 

The treatment for the NK - 2 plane is shown in Figure 12d. The 
steps involved are identical to those described in Figure 12c, except 
that since this is the last value of k for which a plane of predicted 
values is to be calculated no partial-predictor values are computed. 

The final step of the sweeping procedure is shown in Figure 12e, 
First a boundary condition is applied to the predicted values at NK 
creating enough information to compute the completor-corrector values 
at NK - 1. The final plane of corrected values is then obtained from 
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the same boundary condition which was applied to the predicted values. 
Then both n+1 level planes at NK - 1 and NK are stored in the three- 
dimensional solution array by overwriting the final two planes of the 
n level solution. The n+1 level solution is now completely stored in 
the three-dimensional solution array. This completes one sweep of the 
sweeping procedure as well as a single time step. 

The final solution is obtained in a manner typical of the time- 
c endent approach. Using the previously outlined procedure, the 
computation is advanced in time from the initial conditions until the 
solution is changing by a sufficiently small amount. At this point 
the solution is said to have reached ’’steady state." 

C. Shock Jump Conditions 

The bow shock wave is a boundary of the computational domain 
across which the exact shock jump conditions are applied (Rankine- 
Hugoniot equations). Exact treatment of a shock wave by such a pro- 
cedure is called "shock fitting." Because the bow shock position is 
not known a priori, a procedure for determining its position is in- 
corporated with the application of the shock jump conditions. The 
three-dimensional shock fitting procedure used in the present study 
was developed from the two-dimensional procedure used by Tannehill 
et al. [54]. It is similar to the method used by Daywitt and Anderson 
[47] (two-dimensional, time-marching problems) and the methods used 
by Thomas et al. [62] and Sutler et al. [63] (three-dimensional, steady, 
space-marching problems) . 
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The three-dimensional shock jump conditions used in the present 
study are applied only one plane at a time in order to be compatible 
with the sweeping procedure described in the last section. The solu- 
tion variables at i - 1 and along j (l^j ^NJ) for all completor-pre- 
dictor and completor-corrector planes are computed from the shock jump 
conditions. Shock jump conditions as well as other boundary conditions 
are not required for the partial -predictor or the partial-corrector 
steps. A detailed description of the three-dimensional shock jump 
procedure for a perfect gas is now presented. A similar presentation 
involving a real gas in equilibrium is presented in Appendix A. 

The three-dimensional shock jump procedure starts by computing the 
n+1 values for the bow shock radius from the Euler predictor equation 
given by 


r 


n+T 




where r is the shock radius, r~ the radial shock velocity and At the 
s b t 

time increment. Equation 3.9 is used to compute the predicted values 
for the shock radius everywhere except on the boundaries where the 
following boundary conditions are used: 
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The next step involves the calculation of the pressure immediately 

behind the bow shock (p^ . , = p„) . This is accomplished by using the 

partial-predictor, completor-predictor combination given in Equations 

3.3 and 3,4, The resulting predicted values are decoded into the 

primitive variables ( p, u * u Qj u j e) by using Equations 3.7 and 3.8. 

3T y s 

The pressure is then obtained from the equation of state (Equation 2.11). 

Because Equation 3.3 is applied on a boundary (i=l) the backward 
differences used to approximate some of the derivatives in the shear 
stress and heat flux terms require unobtainable information at i - 1. 

To correct this situation forward differences are used wherever back- 
ward differences are not valid. The use of forward differences in 
place of backward differences reduces the method's second-order accuracy 
in the vicinity of the bow shock but does not harm the overall solution 
accuracy. 

The shock slopes (r s and r s ) at the n+1 level are calculated 

H z 

next by using the following central difference formulas 



3.11 


3.12 


where A9 and Az are the physical grid increments in the 0- and z-direc- 
tions, respectively. Values for the shock slopes on the boundaries 
which cannot be computed from either Equation 3.11 or 3.12 are obtained 


by 
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The primitive flow field variables obtained from the modified 
MacCormack 1 s scheme applied at i = 1 are now recomputed from the pressure 
(P 2 )> the shock slopes ( r sg an< ^ r s z ^ arlt ^ the s h° c k jump conditions. 

These conditions are derived from the Rankin e -Hu goniot equations and 
shock geometry considerations in Appendix B and are given by 


P2 = P => 
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3.17 



3.18 


3.19 


where the 2 subscript indicates conditions just downstream of the bow 

shock (l,j,k) at the n+1 level. is the component of freestream 

velocity normal to the bow shock and measured in a coordinate system 

fixed with respect to the boxi; shock (positive inward). All r^, r , 

r c and r Q values are at the (j,k) position in space and the n+1 
t 

level in time. The quantities p , p , y 3 u^ > a nd u (infinity 

n ‘eo’ r <s r,eo u 0, co Zjco 

conditions) are either freestream conditions or flow conditions doxm- 
stream of the planar impinging shock. 

Before Equations 3.14-3.19 can be applied, a test must be made 
to determine which set of infinity conditions are to be used. This 
test must be performed for each grid point at which shock jump condi- 
tions are applied. First, the position of the intersection line be- 
tween the planar impinging shock and bow shock waves must be found. 
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The impinging shock is a fixed plane in the physical domain, and is 
expressed analytically as follows 


r 

P 


j»k 


X int 

z int 


Az(k - 1) + x 


int 


1 

cos[A9(j - 3/2) ] 


3.20 


where r^ is the planar impinging shock radius and x^ t and z ^ n(: are the 
r- and z-intercepts, respectively, of the planar impinging shock in. 
the 0= 0° plane. Of course, the bow shock cannot be described analyti- 
cally. Therefore, the description given by the two-dimensional array 

of values, r e . , , must be used. For a given value of j the intersection 
s J >k 

line position is found by finding the value of k for which 



r 1 < 0.0 

S j,k+U 
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After finding which side of the intersection line that a particular 
bow shock point lies, the correct set of infinity conditions required 
for the shock jump relations are determined. The bow shock moves as 
part of the time-dependent solution. Therefore, the position of the 
intersection line is recalculated before each time step. This com- 
pletes the shock jump predictor step. 

The shock jump corrector step is very similar to the predictor 
step with only two major differences. First, instead of using Equation 
3.9 to compute the corrected shock radius, a modified Euler formula is 


used and is given by 
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Second, when calculating the pressure behind the bow shock. 
Equation 3.5 and a modified version of Equation 3.6 are used which is 
given by 
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n+1 

j> k 
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3.23 


Note that a forward difference is used to approximate the 3F/A§ deriv- 
ative instead of the usual backward difference. This reduces the 
method's second-order accuracy in the vicinity of the bow shock but 
does not harm the overall solution accuracy. 


D . Boundary Cond i ti on s 

The boundary conditions must be incorporated with the sweeping 
procedure presented earlier. Since only one k-plane of predicted or 
corrected values is computed at a time, only one k-p!ane of boundary 
conditions can be completed at a time. Figure 13 shows a schematic of 
a typical lc-plane in the physical domain with the various boundary 
conditions indicated. 

Along the bow shock the shock jump conditions are applied (as 
previously discussed) and are given by Equations 3.14-3.19. The 
wall boundary conditions are determined by specifying an isothermal 
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4 - OUTFLOW BOUNDARY: 
NJ 



Figure 13. Physical 2 -plane showing boundary conditions 
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wall, a zero normal pressure gradient, and the no-slip condition. These 
conditions are given by 


(u ) . . = 0.0 

r NI,j ,k 

<U 0 ) NI,j,k “ °’° 

( VNl,j,k = °-° 


e„„ , , = c T = e = constant 
NI,j,k v w w 

P NI,j,k = P NI-l,j,k e KI-l,j,k/ e w 


3.24 


where T is the specified wall temperature. A perfect gas assumption 
is required to formulate these boundary conditions. 

The 0-outflow boundary is treated with a second-order extrapola- 
tion and is given by 


U NJ ‘ 3 %J- 


u 


NJ 


+ U 


NJ-3 


3.25 


where U stands for any of the primitive solution variables. This 
boundary condition is stable provided the outflow Mach number in the 
inviscid region of the shock layer is supersonic. When this boundary 
condition is used near the intersection line for the shock impingement 
case, numerical difficulties are experienced. The problem arises be- 
cause Equation 3.25 attempts to extrapolate across the nearly discon- 
tinuous region just downstream of the intersection line at i - 1 and 2. 
To correct this a zeroth-order extrapolation is used instead of the 
second-order extrapolation along i = l and 2. 
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A plane of symmetry lies equidistant between the two grid planes 
at j = 1 and 2. This permits a reflective set of boundary conditions 
at j = 1 which is given by 

( u r)i,l,k ” ( u r)i,2,k 

( u e)i,l 5 k = "( u e)i, 2 ,k 

( U z)i,l,k ( U z)i,2,k 

p i,l,k 
e i , l,k 

Notice that all solution components are reflected in a positive manner, 

except the tangential velocity component (U ), which is reflected in a 

0 

negative manner. 

Boundary conditions for the inflow and outflow planes in the z- 
direction are now discussed (see Figure 14). The flow variables at 
the inflow plane (i>jjl) are held fixed for all time equal to the con- 
ditions from a swept infinite cylinder solution calculated prior to the 
shock impingement solution. The flow variables at the outflow boundary 
are determined using a zeroth-ordet extrapolation given by 

U i,j,NK = U i,j,NK~l 3,28 

where U stands for any of the five primitive solution variables. 


P ij2,k 


= e 


i>2»k 


3.27 



IMPINGING 



Figure 14. Physical 9- -»lane showing boundary conditions 
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E. Initial Conditions 


A swept infinite cylinder solution is calculated prior to the 
shock impingement solution with freestream conditions matching the 
conditions behind the impinging shock (region 1, Figure 14). The 
conditions behind the impinging shock are computed using oblique shock 
equations [64 \ 

The infinite cylinder solution is computed from the same computer 
code which computes the three-dimensional shock impingement solutions. 
An infinite cylinder solution, by definition, must have all z-gradients 
equal to zero. Rather than impose this condition directly, by setting 
all z-derivatives equal to zero and thereby changing the computer code 
drastically, it is done indirectly by changing certain boundary condi- 
tions. That is, the z-derivatives are still calculated but automati- 
cally come out to be zero because all z-planes are forced to be identi- 
cal. This identity can be achieved simply by changing two boundary 
conditions. First, reference to the impinging shock is removed at the 
bow shock thus taking away the z-gradients across the intersection 
line. Second, the inflow z-plane (k = 1) , instead of being held fixed, 
is determined by a zeroth -order extrapolation given by 


U. . = U. . 0 


3.29 


where U stands for any of the five primitive solution variables. This 
removes the z-gradient which could exist between the first two z-planes 
of the solution. In addition, all z-planes of the infinite cylinder 
initial condition solution must be identical. With all of these 
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conditions satisfied, a mechanism for producing a nonzero z-gradient 
does not exist. This procedure yields an identical infinite cylinder 
solution in each z-plane of the finite-difference grid. 

Because all z-planes of the infinite cylinder solution are identi- 
cal only a single z-plane is needed to completely specify the solution. 
However, because of the three-dimensional computer code's structure, a 
minimum of four z-planes must be used to compute an infinite cylinder 
solution. Two of these four z-planes are boundary condition planes 
(k = 1 and k = 4) . 

At the beginning of a three-dimensional impingement solution cal- 
culation, the flow field solution in each z-plane is set equal to the 
infinite cylinder solution. The shock radius array (r g ) , the shock 
velocity array (r s ^) , and the shock slope arrays ( r sp an< ^ r s 2 ^ as 
as the flow field variables are initialized in this way. The initial 
condition infinite cylinder solution is also the solution maintained 
as the fixed boundary condition at k-1. 

The procedure used to establish the infinite cylinder initial 
condition solution is similar to the procedure used by Tannehill et al. 
[54]. It involves the approximate curve fit of Billig [65] to get 
the shock shape and slope, the shock jump conditions with r g set equal 
to zero to get the flow variables behind the bow shock, a Newtonian 
pressure distribution and wall boundary conditions to get the flow 
variables on the body, and a linear variation between the bow shock 
and the body to get the interior flow field. 
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F. Stability 

To insure numerical stability the time step associated with 
MacCormack 1 s method is limited in size by the CFL condition [61,66], 
which in three-dimensions is given by 


At £ C 


3,30 
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where a is the speed of sound and C is an adjustable constant* 

Numerical solutions involving steep flow field gradients have used 
special artificial smoothing terms to remain stable [67-70]. A modi- 
fied version of the fourth-order product-smoothing term first developed 
by MacCormack [71] is used in the present study. The original smooth- 
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where the U’s are defined by Equation 2.27, the summations on l indi- 
cate one term for each of the three spacial directions (i,j,k), the 
c's are constant coefficients, the K's are variable coefficients de- 
fined by 
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n = l p i+2,_i,k ~ 2p i+l, i,k + P i,i,k _ 
( p i+2, j ,k + 2p i+l s j,k P i»j,k/ 


K M-1 
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and the 6 () and A.Q operators are forward and backward differences* 

Jo XJ 

respectively, defined by 
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6 i U “ U i+l,j,k U i s j jk 
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When smoothing is desired, S? . . (Equation 3.31) is added to the pre- 

J 3 5 K 

dictor step (Equation 3.1), and S. 11 *, 1 . (Equation 3.32) is added to 

1 3 3 

the corrector step (Equation 3.2). 

The variable K-coefficients are composed of a normalized second- 
order difference of flow field density. The K subscript indicates 
which direction the difference is in (i,j or k) and also the center 
of the difference. The K superscript indicates the time level used 
for the values of density. This always positive variable coefficient 
is essentially zero in the smooth regions of the flow field and ap- 
proaches a maximum value of one in regions of large point-wise oscil- 
lations. The theoretical maximum value of the coefficient product in 
front of the U-dif ferences, namely the cK products, is one-half. The 
K-coefficients can theoretically reach a value of one, causing the c^ 

c , and c, constant coefficients to be restricted to a value of one- 
k 

half or less. However, since the K-coefficients in practice are 
usually much smaller than one, the constant coefficients can be much 
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larger than one-half. Of course, smaller amounts of smoothing or 
smoothing in only one or two of the three spacial directions can also 
be obtained. This is achieved by independently setting the constant 
coefficients in each of the three directions equal to the appropriate 
values, 

A major difficulty arises when applying this smoothing method to 
the sweeping procedure presented earlier. To obtain both the predicted 
and corrected smoothing terms at a particular value of k, solution vari- 
ables from five k-planes are needed (k + 2, k+ 1, k, k - 1, k - 2) . Pro- 
viding the information from such a large number of k-planes does not 
fit within the framework of the sweeping procedure. At most, infor- 
mation from only three k-planes can be supplied. This leads to a 
modification of MacCormack' s smoothing technique which is given by 


S? P . V 


?. c K 


A=i, J 


mV 


n 


W 


n) 


4- c 


s cp 

i) 



3.36 


3.37 


S 


c 

j jk 


1 

2 



K 


i - 1 


¥ 







3.38 


where all notation is as it was previously defined. “When smoothing is 
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desired, S?^. , , the partial-predictor-smoother step, is added to the 
1 ) J > “ 

cp 

partial-predictor step (Equation 3.3), S. . , , the completor-predictor- 

1 > J > k 

smoother step, is added to the camp letor -predictor step (Equation 3.4), 

Q 

and S. . . , the corrector-smoother step, is added to the completor- 

1 1 i k 

corrector step (Equation 3.6). 

Note that the modification does not alter the terms in either the 
i- or j-directions, only the k-direction. Instead of using K-coeffi- 
cients centered at three different k-locations (k+1, k, k - 1) , all IC's 
in the k-direction are centered at k. This modification alone reduces 
the number of k-planes involved to three. 

All quantities in the partial-predictor-smoother step (Equation 

3.36) are obtained from two k-planes (k and k - 1) , except p” . , 

1 j J > k+ I 

which appears in the coefficient. Since only the k and k - 1 planes 
of the n-level solution exist in SCM at this time, extra density in- 
formation at k+1 must be obtained from LCM. These coefficients are 
then saved and used in the completor-predictor-smoother step (Equation 

3.37) . 

Instead of breaking the corrector -smoother step into too parts, 
it ii applied as one step just after the completor-corrector step. 

This is because a third plane of densities (at k + I) cannot be read 
into SCM from LCM to compute a partial-corrector-smoother step. Because 
of the sweeping procedure, the n + 1 - level solution at k+1 has not 
been computed yet. At the time the partial-corrector step is completed, 
the entire n + 1 solution at k-1 is saved in SCM. This can be con- 
sidered as a partial-corrector-smoother step even though no computations 
are actually made. Then, after Lhe completor-corrector step is computed, 
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the entire corrector-smoother step (Equation 3.38) is computed with all 
three k-p lanes existing in SCM. 

When Equation 3.12 is used to compute the shock slopes in the z- 

direction (r s ) an unstable point-wise oscillation results in the z- 
z 

direction. This instability is strongest near the z-direction inflow 
plane. To alleviate this problem a ohock radius smoother is applied 
during the corrector step which is given by 
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where a is an adjustable constant. When a has a value of 0.25 the 

amount of smoothing which occurs corresponds to the optimal flattening 

of a sawtooth curve (r va z distribution) into a straight line. Larger 

s 

values of a will produce "over smoothing" and if large enough (above 
0.5) will actually feed the instability. Values of a smaller than 0.25 
(up to an order of magnitude smaller) will allow oscillations to exist 
but will still keep the procedure stable. 
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IV. DISCUSSION OF RESULTS 

The present method was used to compute two-dimensional blunt body, 
swept infinite cylinder, and three-dimensional shock impingement solu- 
tions, Results for these three problems are presented in this chapter 
along with the computational statistics of the computer code. In all 
cases considered in this chapter, air is assumed to be a perfect gas.. 

The Prandtl number (Pr) , the coefficient of specific heat at constant 

pressure (c ) , and the ratio of specific heats (Y) are held constant and 
P 

are given by 

Pr = 0.72 

c = 1004.58 m 2 /sec 2 - °K 4.1 

P 

y = 1.40 

The simplest case, that of the two-dimensional blunt body solution, is 
presented first. 


A. Two -Dimensional Blunt Body Solution 


The two-dimensional blunt body solution was computed using the three- 
dimensional code in the same manner as the infinite cylinder solutions 
but with the freestream cross-flow velocity set equal to zero. 

The freestream conditions chosen for this two -dimensional computa- 
tion were 


M - 5.94 



p = 545.8 N/m 2 

CD 

T = 62.8 °K 


= 17326 


4.2 
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where M is the Mach number, p is the pressure, T is the temperature, and 
Re^ is the Reynolds number based upon the cylinder diameter. Other con- 
ditions pertinent to this case were chosen to be 

D = 0.025 m T = 1444.4 °K 

W 4.3 

X = 0 ° p « 0 

where D is the cylinder diameter, X is the cylinder sweep angle, 8 is 

the radial mesh stretching parameter, and T is the constant wal;. tern- 

w 

perature. Note that a X of aero indicates no sweep and therefore, a 

two-dimensional blunt body case, A 8 of sero indicates no stretching 

of the mesh in the radial direction. The adiabatic wall temperature 

(T ) in the stagnation region was computed to be 504.9 °K from 
3W 

T a W - T e (l+r 

where e denotes the edge of the boundary layer and r is the recovery 

factor set equal to 0.85 for a laminar boundary layer. On comparing the 

T and T values, it is apparent that this case has an extremely hot 
w aw J 

wall. 

The computational mesh consisted of 21 equally spaced points in 
both the radial and tangential directions. The tangential outflow 
boundary was located along a radial ray 83.8 ° above the stagnation 
streamline. The results of this computation are shown in Figures 15 and 
16. 

A plot of the pressure distribution around the cylinder is shown in 
Figure 15. The circles are the current numerical results. The solid 
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Figure IS. Comparison of pressure distributions on a two-dimensional 
circular cylinder. 



Figure 16. Comparison of shock shapes around a two-dimensional circular 
cylinder. 


DISTANCE, r , cm 
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curve is a fairing of experimental data due to Beckwith and Cohen [72]. 
This faired experimental curve consists of data from both two-dimensional 
blunt body and infinite cylinder results over a range of Mach number. 

This pressure distribution is independent of sweep angle and normal Mach 
numbers above two. The dashed curve shown in Figure IS is the result 
predicted by the modified Newtonian pressure distribution given by 


■ - (l - cos^0 + ^ 4.5 

p V p P / 

r stag r stag r stag ' 


where P sta g is the stagnation point pressure computed from the Rayleigh 
pitot formula [64], 

All three results are in excellent agreement in the stagnation re- 
gion of the blunt body. In fact, the error between the present numeri- 


cal result and the stagnation value of pressure as predicted by the 
Rayleigh pitot formula is less than 0.1 percent. As expected the modi- 
fied Newtonian pressure disagrees with the other two results away from 


the stagnation region. The faired experimental curve and the numerical 
results start disagreeing at about 9 = 50 °. This disagreement in- 
creases to a maximum of about five percent error at 0 = 80 °. The 


scatter of the experimental data used to obtain the faired experimental 
curve in Figure 15 was about ± five percent at R = 80 °, 


A comparison of the numerically predicted shock shape with the em- 
pirical result of Billig [65] is shown in Figure 16. The circular sym- 
bols represent the numerical result while the empirical result is shown 
as a solid curve. The error between the two results, to a large extent, 
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is due to the extremely hot wall. The large wall temperature heats the 
boundary layer causing the fluid density near the wall to decrease. 

This apparently has the effect of increasing the bow shock standoff 
distance as seen in the next section on infinite cylinder solutions. 

B. Swept Infinite Cylinder Solutions 

The swept infinite cylinder solution (as already mentioned) was com" 
puted using the three-dimensional computer code by restricting all s-gra- 
dients to be zero. This was accomplished with the appropriate choice of 
boundary and initial conditions. Three swept infinite cylinder solutions 
were computed to establish the validity of the three-dimensional computer 
code and are presented in this section. In addition, the swept infinite 
cylinder solution used as the initial condition of the three-dimensional 
shock impingement case is also presented. In all of these cases, the 
computational mesh consisted of 21 mesh points in both the radial and 
tangential directions. The tangential outflow boundary was in all cases 
located along a radial ray 83.8° above the stagnation streamline. 

The specified conditions chosen for the three infinite cylinder 
solutions are given in Table 1. The freestream Mach number has been 
divided into two components, the normal component (M^ m ) and the cross- 
flow component (M ) . Except for the addition of a cross-flow compo- 

z , co 

nent of velocity, the case A conditions are identical to the two-dimen- 
sional case. Because of this, the previously computed two-dimensional 
solution was chosen a ; the initial condition for the case A infinite 

cylinder solution. The time-dependent development of the u -velocity 
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Table 1. 

Specified conditions 
solutions. 

for the infinite cylinder 

Specified 

Quantity 

Case A 

Case B 

Case C 

M 

CO 


6.45 

6.45 

6.45 

M 

r,« 


5.94 

5.94 

5.94 

M 


2.51 

2.51 

2.51 

Re 

D, co 


18816 

18810 

18816 

p (U/m 2 ) 

CO 


545.8 

545.8 

545.8 

T (°K) 

CO 


62.8 

62.8 

62.8 

X (rteg) 


25.0 

25.0 

25.0 

D (m) 


0.025 

0.025 

0.025 

T (°K) 


1444.4 

1444.4 

411.8 

T 

aw, stag 

<°K) 

573.2 

572.0 

572.5 

P 


0.0 

1.2 

1.2 
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profile along the stagnation line for case A is shown in Figure 17, 

The u -velocity profile develops from the bow shock toward the body. 
z 

Moderate oscillations occur in the radial direction during the time- 
dependent process but completely disappear as the steady-state solution 
is approached. Initially, the solution moves rapidly toward sceady 
state, however, the final small oscillations damp out very slowly as 
steady state is approached. The process essentially does not disturb 
the initial solution in the r0-plane. This suggests independence be- 
tween the two- dimensional solution and the cross-flow direction solution. 
The second infinite cylinder solution (case B) was computed with 
exactly the same conditions as the first, but with a different radial 
distribution of mesh points. The mesh was refined near the body with 
the stretching parameter (3) equal to 1.2. For this value of 0, the 
first grid point off the body is located at 2.3 percent of the shock 
standoff distance as compared with five percent for an equally spaced 
mesh. The mesh distribution for 21 points is tabulated in Table 2 for 
several values of p . The distance between the body and the bow shock 
has been normalized to one. The first and second derivatives of the 


stretching transformation are also tabulated in Table 2. 

The u^ -velocity and temperature profiles near the stagnation region 

(0 = 0°) and the u -velocity profile at f) = 36.5° are shown in Figures 

0 

18, 19 and 20, respectively. The distance between the bow shock and 
the body (x) is normalized to one for all profiles. For comparative 
purposes the solutions for no stretching (case A) and for stretching 


(case B) are both plotted. All three profile comparisons are in excel- 


lent agreement. 
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Figure 17. Time development of the u z -velooity profile near the stagna- 
tion region of a swept infinite cylinder. 
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Figure 18. u s -velocity profile near the stagnation region of a swept 
infinite cylinder. 



NORMAL DISTANCE TO BODY 


76 



TEMPERATURE, T, °K 


Figure 19. Temperature profile near the stagnation region of a swept 
infinite cylinder. 
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Figure . 20. ug-velocity profile at 0 - 36.5 on a swept infinite 
circular cylinder. 
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The u -velocity profile (Figure 18) is typical of all such profiles 

on the infinite cylinder regardless of the tangential position. The 

inviscid portion of the profile has a constant value equal to the a- 

component of freestream velocity. The boundary layer portion of the 

profile changes only slightly with tangential position. This occurs 

in this computation because the mesh expands (due to the shock layer 

growth) at the approximate growth rate of the boundary layer. 

The temperature profile (Figure 19) shows the extent of the hot 

wall. The u -velocity profile at 0 = 36.5 ° (Figure 20) has a velocity 

0 

overshoot which is characteristic of such a hot wall solution. The 
better boundary layer resolution of the stretched solution is apparent 
in all three of the profile comparisons. However, even better resolu- 
tion through more grid points or larger amounts of stretching would be 
desirable fcr this solution. 

The third infinite cylinder solution (case. C) was computed with 


exactly the same conditions as case B but with a wall temperature 
slightly less than the adiabatic wall temperature (see Table 1). There- 
fore, case C is a cold wall case. As pointed out earlier, the hot wall 
condition seemed to cause a slight increase in the shock standoff dis- 
tance. A plot of the three infinite cylinder shock standoff results 
along with the empirical result of Billig [65] is shown in Figure 21. 

The empirical result was computed from ths normal component of the 
freestream Mach number and therefore, is identical to the two-dimensional 
shock standoff curve. If a hot wall increases the shock standoff dis- 
tance, then a cold wall should decrease the shock standoff distance. 
However, the hot wall cases considered in the previous solutions were 
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Figure 21-. Comparison of shock shapes around a swept infinite 
circular cylinder. 
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extreme cases. The cold wall case under current consideration is only 
a slightly cold wall being very close to the adiabatic wall temperature. 
Therefore, the close agreement shown in Figure 21 between the cold wall 
shock shape and the empirical curve fit is expected. 

A plot of the pressure distribution around the swept infinite cylin- 
der is shown in Figure 22. The three numerical solutions are represented 
by symbols. The solid curve is a fairing of experimental data due to 
Beckwith and Cohen [7 2 1 which was discussed earlier. The agreement is 
generally quite good over the entire curve. 

A u -velocity profile comparison is shown in Figure 23. The hot 
0 

wall profile from Figure 20 (case B) and the cold wall profile of case C 
are both plotted. Both profiles are from the same position on the 
cylinder ( 0 = 36.5°) and have identically stretched meshes. As expected, 
the velocity overshoot has disappeared from the cold wall profile. Other- 
wise, the two curves are almost identical. 

The infinite cylinder solution used as the initial condition for 
the three-dimensional shock impingement solution is presented next. 

The freestream conditions chosen for this case were 
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Other conditions pertinent to this case were chosen to be 
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Figure 22 


Comparison of the pressure distributions on a swept infinite 
circular cylinder. 
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Figure 23. u Q -velocity profile at 0—36.5 on a. swept infinite cylinder 
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These conditions were selected to permit a comparison with the experi- 
ment of Keyes and Hains [2]. The inflow plane in the z-direction for 
the three-dimensional shock impingement solution is fixed at the initial 
infinite cylinder solution. Therefore, the freestream conditions for 
the infinite cylinder solution (except Re„ ) were computed to agree 
with the experimental conditions on the downstream side of the impinging 
shock. These conditions were computed from the given experimental free- 
stream conditions and the oblique shock relations T64]. The freestream 
viscosity was chosen to be an order of magnitude larger than in the 
experiment, thus making the Reynolds number ten times smaller. This 
was done to physically thicken the boundary layer and make its resolu- 
tion possible with fewer grid points. The results of this computation 
are shown in Figures 24-30. 

A plot of the maximum and the minimum radial shock velocity (r St > 
as a function of time step number (n) is shown in Figure 24. These 
velocities should approach zero as the steady-state solution is ap- 
proached. The solution moves rapidly toward a steady state initially 
with large time-dependent oscillations appearing. In the final part of 
the solution, the convergence is very slow and has a monotonic behavior. 

Plots of the numerical wall pressure and the heat transfer distri- 
butions around the swept infinite circular cylinder are shown in Figures 
25 and 26, respectively. The pressure distribution in Figure 25 is com- 
pared with the faired experimental pressure distribution due to Beckwith 
and Cohen [72]. The comparison is in general quite good. No experimen- 
tal results for the heat transfer distribution are available for the 
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Figure 24 . Radial shock velocity as a function of time step number 
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Figure 26, Heat transfer distribution on a swept infinite circular 
cylinder. 
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conditions of the current solution. However, the stagnation line value 
of heat transfer can be compared with theoretical results. 

The stagnation line heat transfer for a laminar boundary layer on 
a yawed infinite circular cylinder can be obtained from [73] 



where the subscript :: stag" stands for values measured along the stagna- 
tion line, k is the free stream coefficient of thermal conductivity, 

CO 

is the coefficient of viscosity computed from the Sutherland viscosity 
law, and T is given as before from Equation 4.4, The value of stagna- 
tion line heat transfer computed from this formula is about 15 percent 
lower than the corresponding numerical value. This may be the result 
of not obtaining enough resolution in the boundary layer. The numeri- 
cal values of wall heat cransfer are computed by evaluating the slopes 
of the total enthalpy profiles at the wall. The slopes are approximated 
by three values of total enthalpy (i = NI, NI - 1, and NI - 2) , This ap- 
proximation is inaccurate for coarse mesh spacings. 

Another consideration for the difference between the two results is 
the accuracy of the theoretical result. Beckwith and Gallagher [73] com- 
pare this theory with experimental results but for Reynolds numbers two 
orders of magnitude larger than the present value. In addition, only 
a limited number of comparisons were made between theory and experiment. 
In light of these . nsi derations, the comparison between the theory and 
the present numerical result seems reasonable. 
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Contour plots of Mach number, density, pressure, and temperature 
for the present infinite cylinder solution are shown in figures 27-30, 
respectively. These contour plots have been drawn by a computer plotter 
and are displayed in the physical plane. No characteristic blunt body 
sonic line is shown in Figure 27 because the z-component of Mach number 
in the shock layer is by itself supersonic. The boundary layer can 
clearly be seen in Figures 27, 28 and 30. Note the unusually large 
boundary layer thickness in the stagnation line region which is charac- 
teristic of an infinite cylinder solution. The general validity of the 
zero normal pressure gradient assumption can be seen from the pressure 
contours shown in Figure 29. The stagnation line region of Figure 29, 
however, does seem to suggest a nonzero normal gradient in pressure. 

This concludes the infinite cylinder solution presentation. The 
results for the three-dimensional shock impingement solution are pre- 
sented in the next section, 

C. Three-Dimensional fhock Impingement Solution 

An infinite cylinder solution was used to obtain the initial con- 
ditions for the three-dimensional wing -leading-edge shock impingement 
solution. The flow field variables in each z-plane of the three-dimen- 
sional initial solution were set equal to the flow field variables from 
the infinite cylinder solution. The planar impinging shock was intro- 
duced by discontinuously changing the freestream conditions across the 
intersection line. The solution was advanced in time by using the 
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numerical procedure previously described until the wall pressures were 

changing by sufficiently small amounts. 

The freestream conditions chosen for this three-dimensional shock 

impingement computation were 

M -5.38 p- 559.1 N/m 2 

06 w 

M *2.51 T * 59.6 °K 4 * 9 

2 , » ® 

Re_ = 19344 0 = 10 ° 

D,» 

where 0 is the angle that the planar impinging shock makes with respect 
to the freestream velocity vector. For this set of freestream condi- 
tions, the pressure ratio across the impinging shock (PR) has a value 
of 3,6. The flow field conditions, which result after the freestream 
flow passes through the impinging shock, must be identical to the free- 
stream conditions used for the infinite cylinder solution. Other condi 
tions pertinent to this case were 
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X = 25° 
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T = 394 °K 
w 

T * 460 °K 
aw 

a = 0.25 

c i = °j = C k * 1,0 


4.10 


where a is the smoothing coefficient for the (previously presented) bow 
shock smoother (see Equation 3.39) and c^, c^, and c^ are the coeffi- 
cients for the 4th-order product smoothing terms used in conjunction 
with the (previously presented) finite difference method (see Equations 
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3.31 and 3.32). Hie sweep angle (A) given in Equation 4.10 corresponds 
to a Type V interaction. 

The three-dimensional computation was broken up into two parts. 
First, a coarse solution was computed in which the computational mesh 
consisted of 21 grid points in all three of the spacial directions. A 
plot of wall pressure at two different stagnation line locations (mid- 
point and outflow) versus time step number for this coarse solution is 
shown in Figure 31. Note the large fluctuations which occur initially 
in the solution. The solution appears converged after approximately 
400 time step iterations. 

Second, a refined solution was computed in which the computational 
mesh consisted of 21 points in both the radial and tangential directions 
and 41 points in the cross-flow direction. The coarse solution was used 
to obtain an initial condition for the refined solution. A linear in- 
terpolation was used to generate an additional 2 -plane of values between 
each already existing z-plane. Thus, the size of the physical domain 
in going from the coarse solution to the refined solution did not change. 
Instead, the number of grid points in the z-direction was doubled (A z 
was halved) to improve the solution. No other changes were made. A 
plot of wall pressure at two different stagnation line locations (mid- 
point and outflow) versus time step number is shown in Figure 32. The 
large fluctuations in the initial portion of the solution are much 

j 

smaller than for the previous case with the coarse mesh. This is 
probably due to the more exact initial condition solution. The solu- 
tion appears converged after approximately 400 iterations. 
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Figure 31 . Time -dependent variation of wall pressure for the 21 x 21 x 21 mesh. 




Figure 32. Time -dependent variation of wall pressure for the 21 x 21 x 41 mesh. 
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Because of the bow shock smoothing technique employed in the three- 
dimensional case, the radial shock velocities (r St > were not able to 
approach zero in the vicinity of the intersection line. The shock 
smoothing technique, which removed erroneous numerical oscillations 
from the bow shock in the z-direction, also removed all desired physi- 
cal "kinks,” The physics of the three-dimensional problem tried to 
restore the physical kinks by developing large shock velocities. The 
bow shock, of course, did approach a steady-state position as the non- 
zero shock velocities and the shock smoothing approached equilibrium. 

The conditions given by Equations 4.9 and 4.10 were all chosen 

(except for Re_ ) to agree with the experiment of Keyes and Hains [2], 

D , « 

As with the infinite cylinder solution, the viscosity was chosen to be 
an order of magnitude larger to physically thicken the boundary layer 
and make its resolution possible with fewer grid points. In addition 
to increasing the boundary layer and shear layer thicknesses, this 
change also affected the wall shear stress and heat transfer values. 
Therefore, nondimensional ratios should be used when comparing these 
quantities with the experiment. 

A comparison of the stagnation plane shock shapes is shown in 
Figure 33. The solid curve is the experimental result of Keyes and 
Hains [2] taken from a Schlieren photograph. The symbols represent the 
present numerical results and are plotted in such a manner that the 
experimental and numerical positions of the impinging shocks agree 
exactly. The experimental results were obtained by allowing a planar 
impinging shock to strike the shock layer on a finite swept cylinder. 

The intersection point along the stagnation plane was only three centi- 
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meters downstream from one end of the cylinder. The shock standoff 
distance for the initial numerical z-plane is therefore, much different 
than the corresponding value of the experimental result. When these 
curves are examined in light of this difference the comparison seems 
quite good. Note the effect of the bow shock smoothing in the numeri- 
cal result. The physical kinks which appear in the experimental Schlleren 
photograph did not develop to the full extent in the numerical result, 

A comparison of the stagnation line wall pressure is show in 
Figure 34. The symbols represent the experimental results of Keyes and 
Hains [2] and were obtained from static pressure ports in conjunction 
with electrical strain-gage pressure transducers. The solid curve 
represents the present numerical result. The z-direction position of 
the numerical pressure distribution curve was determined from the shock 
shape comparison of Figure 33. The general trend of the comparison is 
quite good. However, the peak value in the experimental curve, which 
is caused by a boundary layer interaction with the transmitted shock, 
is not reproduced in the numerical results. A small peak does occur in 
the numerical results indicating that the transmitted shock is partially 
formed. 

A comparison j f the stagnation line heat transfer is presented in 
Figure 35, The symbols represent the experimental results of Keyes 
and Hains [2] obtained from the phase -change coating technique. The 
solid curve is the present numerical result. A peak in the heating rate 
is measured for both the numerical and experimental results although 
the heights of the peaks are not in good agreement. The coarse grid, 



PRESSURE RATIO 



100 




HEAT TRANSFER RATIO, qyq 


2.4 



Figure 35. Stagnation line heat transfer rates. 
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numerical smoothing, and increased physical viscosity probably all 
contributed to the poor resolution of the transmitted shock and there- 
fore, to the poor agreement. 

More insight into the numerical solution can be obtained from con- 
tour plots of the important solution variables. Contour plots of Mach 
number, density, pressure, and temperature are shown for selected z- 
plane and 0-plane locations in Figures 36-50. Each of these contour 

plots was drawn with a computer plotter. The Mach number, density, 

-2 . 3 

pressure, and temperature increments were 0.05, 1.03x 10 kg/m , 

2 

957.6 W/m , and 11.11 K, respectively. The sketch in the upper left- 
hand corner of each contour plot figure identifies the approximate lo- 
cation of the plane being plotted. For clarity all plots are shown in 
the physical domain. 

Mach number contour plots at k = 8, 14, 20, 28, and 36 are shown 
in Figures 36, 37, 38, 42, and 43, respectively. This series of contour 
plots shows the effect of shock impingement at various z-plane locations. 
The stagnation plane impingement point is between k=ll and 12. There- 
fore, the Mach number contour plot at k = 8 (Figure 36) shows no sign 
of the impinging shock. The next five Mach number contour plots snow 
the development of a strong shear layer which is the result of the 
impinging shock. The shear layer, which is seen in cross-section, em- 
anates from the bow shock and approaches the body. 

In addition to Mach number, z-plane contour plots of density, 
pressure, and temperature at k = 20 are shown in Figures 39-41, respec- 
tively. The shear layer is apparent in both the density and temperature 
contour plots but not pressure. An intersection point between the 
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Figure 36. z-plane Mach number contours for a swept circular cylinder under 
the influence of an impinging shock (k = 8) . 
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Figure 37. z-plane Mach number contours for a swept circular cylinder 
under the influence of an impinging shock (k = 14). 



Figure 38. z-plane Mach number contours for a swept circular cylinder 
under the influence of an impinging shock (k = 20). 
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M =5.94 



Figure 39. z-plane density contours for a swept circular cylinder 
under the influence of an impinging shock (k = 20). 
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Figure 41. z-plane temperature contours for a swept circular cylinder under 
the influence of an impinging shock (k = 20). 
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Figure 42. z-plane Mach number contours for a swept circular cylinder under 
the influence of an impinging shock (k = 28). 
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Figure 45. 0-plane density contours on a swept circular cylinder 
under the influence of an impinging shock (j = 2). 
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Figure 46. 0-plane pressure contours on a swept circular cylinder under the 
influence of an impinging shock ( j = 2) . 
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Figure 47. 0-plane temperature contours on a swept circular cylinder under 
the influence of an impinging shock (j = 2). 
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planar impinging shock and the bow shock i. e quite apparent in several 
of the z-plane contour plots (k=l4, 20, and 28). Of course, no inter- 
section point appears in the z-plane contour plots at k = 8 or 36, be- 
cause the range of the intersection line only extends from k = 12 to 
k = 30. The effect of the zeroth-order extrapolation on the 0-outflow 
boundary condition at i ~ 1 and 2 is seen in most of the z-plane contour 
plots (bow shock vicinity of the 0-outflow) . As pointed out earlier, 
this boundary condition was employed to prevent extrapolation across 
the sharp gradient in the flow near the intersection line. Although 
it succeeds in this respect, it also introduces an inconsistent region 
in the flow field. 

Mach number contour plots at 0 = 2.1° (j = 2, stagnation plane), 
10.7° (j = 4), 27.9° (j=8), and 62.3° ( j = 16) , are shown in Figures 
44, 48, 49, and 50, respectively. The flow field structure is best 
studied from these 0-plane contour plots. The strong shear layer pro- 
duced by the shock impingement appears as a dark band of coalesced con- 
tour lines. This shear layer emanates from the bow shock near the in- 
tersection line, approaches the body, and finally, passes through the 
s-outflow boundary. 

A secondary coalescence of contour lines (Figure 44) which repre- 
sents another (weaker) shear layer, lies between the bow shock and the 
first shear layer. An inspection of the numerical output shows that 
the flow field between these two shear layers is subsonic. This region 
of subsonic flow enists only for small values of 0. The extent of the 
subsonic flow region can be seen in the z-plane Mach number contour 
plots (Figures 37, 38, 42, and 43). Further inspection of the 
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numerical output: reveals a weak jet existing just above the weak shear 
layer. Both the jet and the weak shear layer follow the stronger 
shear layer (in a parallel manner) through the z-outflow boundary. 

In addition to Mach number, 0-plane contour plots of density, 
pressure, and temperature at 0 = 2,1° (j = 2, stagnation plane) are 
shown in Figures 45-47, respectively. As expected, the strongest shear 
layer is apparent in both the density and temperature 0-plane contour 
plots but not the pressure. 

A feature of the flow field which is apparent in the density and 
pressure 0-plane contour plots is the initial formation of the trans- 
mitted shock. The darkest coalescence of density or pressure contours, 
which is approximately perpendicular to the body represents a very 
sharp expansion. The vertical contour lines immediately downstream of 
the sharp expansion are the beginnings of a transmitted shock. Still 
further downstream is another (weaker than the first) expansion. The 
sharp expansion, followed by a slight pressure jump and then another 
expansion can best be visualized by referring back to Figure 34 which 
is a plot of stagnation line wall pressure. 

Other noteworthy aspects of the 0-plane contour plots include the 
boundary layer, the thickening of the shock layer for increasing values 
of 0, and the approximate validity of the zeroth-order extrapolation 
for the z-outflow boundary condition. 

A Schlieren photograph of a Type V shock impingement on a swept 
circular cylinder is shewn in Figure 51. This experimental result is 
due to Keyes and Hains [2]. The conditions given in Equations 4.9 and 
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Figure 51. Schlieren photograph of a type V shock interference (Keyes and Hains) 
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4.10 (except, of course, the freestream Reynolds number) were chosen 
to match this experiment. Note the existence of the transmitted shock, 
the shear layer, and the jet. The development of the shock layer from 
the end of the cylinder is also shown in this photograph. 

The failure of the numerical method to more adequately predict the 
transmitted shock is due largely to the coarse mesh. Another contrib- 
uting factor is the bow shock smoothing which keeps the two bow shock 
kinks from forming. This in turn directly affects the shock layer de- 
tail especially the transmitted shock. 

Contour plots of the wall pressure and heat transfer are shown in 
Figures 52 and 53, respectively. The overall effect of the shock im- 
pingement at the wall is seen quite vividly from these figures. The 
peak in the heat transfer moves around the body in much the same 
fashion as the intersection line on the bow shock. The effects of the 
shock impingement are largest at the stagnation line and diminish in 
the tangential direction until little change from the infinite cylinder 
solution is experienced at the ft-outflow boundary. 

D. Computational Statistics 

The two most important computational statistics associated with 
the solution of any problem by a finite-difference technique are the 
execution time and the storage requirements. The execution time, for 
the most part, determines the cost of the computation. The storage re- 
quirement determines the size of the computer needed for the computa- 
tion. The storage requirement directly limits the computational mesh and 
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very possibly the quality of the solution (especially for three-dimen- 
sional problems). Statistics such as these are examined in this 
section. 

The computational statistics of the solutions presented in this 
report are given in Table 3. The number of time steps per case varies 
over a large range because of the different levels of desired conver- 
gence. For instance, the infinite cylinder solution (used as the ini- 
tial condition for the three-dimensional shock impingement solution) 
was carried to a much higher level of convergence than the other infi- 
nite cylinder cases. 

As previously mentioned, the large core memory (LCM) is used to 
store the three-dimensional solution array. Therefore, the required 
amount of LCM is proportional to the grid size. The small core memory 
(SCM) is used to store the program and the intermediate result arrays. 
The amount of storage listed in Table 3 does not include buffer storage 
or any storage allotted to the system. 

The final statistic shows how much execution time is required on 
a per time step per grid point basis. The last two values in this 
category (three-dimensional solutions) are higher than the others 
because proportionally they contain a smaller number of boundary con- 
dition points. All execution times are for a Control Data Corporation 
7600 Computer. 

The small amount of storage required for the three-dimensional 
solution (only 131K words) , is an achievement since the solution array 
by itself required 90K words. With the savings in storage produced by 



Table 3. Computational statistics 


Two- 

Dimensional 

Infinite Cylinder Cases 


Three -Dimen s Iona 1 


Initial 

Shock Impingement 

Case 


Condition 

Coarse 

Fine 


Case A Case B Case C 

Case 

Mesh 

Mesh 


Time steps 

900 

700 

800 

No. of runs 

6 

6 

3 

Storage 
SCM (K words) 
LCM (K words) 

30 

11 

30 

11 

30 

11 

Mesh 

dimensions 

21x21x5 

21x21x5 

21x21x5 

Execution 
time (sec) 

475 

371 

410 

Execution time 
per time step 
per grid pt. 

0.24 

0.24 

0.23 


(msec) 


600 

1500 

750 

400 

3 

14 

14 

13 

30 

35 

41 

41 

11 

9 

46 

90 

21x21x5 

23x21x4 

21x21x21 

21x21x41 

310 

550 

2832 

2797 

0.23 

0.21 

0.41 

0.39 
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the sweeping procedure, larger more productive computational meshes can 
be used. 
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V. CONCLUDING REMARKS 

A time-marching finite-difference method has been used to compute 
the three-dimensional wing- leading-edge shock impingement problem. 

The bow shock was treated as a sharp discontinuity across which the 
exact shock jump conditions (Rankine-Hugoniot relations) were applied. 
The impinging shock (assumed planar) was introduced by discontinuous ly 
changing the freestream conditions across the intersection line at the 
bow shock. The compressible Navier-Stokes equations were used for this 
three-dimensional computation. The present method does not require 
any prior information about the shock impingement flow field to be com- 
puted , as is the case with previous semiempirical approaches. In 
addition, since the shock layer flow field is automatically "captured" 
in the same manner in each computation, it is possible, in principle, 
to compute all three types of wing- leading-edge shock impingement 
with the same computer code. 

A special storage- saving procedure for sweeping through the finite- 
difference mesh has been developed. This sweeping procedure reduces 
the required amount of computer storage by at least a factor of two 
without sacrificing the execution time. This savings in storage allows 
for larger more productive computational meshes which is especially 
useful for three-dimensional computations. 

The shock impingement solution presented in this study demonstrates 
the feasibility of three-dimensional time-dependent computations in- 
volving the Navier-Stokes equations. At the time of this writing, only 
one other time-dependent solution of the three-dimensional Navier-Stokes 
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equations is known to exist T74], The computational demands associated 
with this problem, however, affect the quality of the present three- 
dimensional solution. The artificially reduced Reynolds number and 
the relatively coarse mesh are both limitations resulting from these 
computational demands. As a result of this study, however, it is felt 
that these limitations can be removed with the present method and with 
currently available advanced computers. 
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VI. RECOMMENDATIONS FOR FURTHER STUDY 

Several recommendations concerning the extension of this work can 
be suggested. First, a new version of the shock smoothing technique 
(the old version is given by Equation 3.39) could greatly improve the 
results. Obviously this new technique should eliminate the oscilla- 
tions in the bow shock, but, should not interfere with any physical 
"kinks" which might appear in the bow shock. 

The alleviation of the bow shock oscillations could also be 
achieved in the shock slope calculation without the addition of any 
numerical smoothing. For instance, the use of a one-sided difference 
formula in place of the central difference formula to compute the shock 
slope might relieve the numerical instability causing the oscillations 
and therefore, remove the need for a bow shock smoother. One-sided 
formulas which could be used for this purpose are given by 

f r s ^ = ( r s " r s ) /Az 6,1 

\ z/j ,k \ j,k+l j,k/ 


(x \ = /4r - r - 3r \/2Az 6.2 

\ S zjj,k \ S j,k+1 S j,kf2 s j,ky 

where Equation 6.1 is a first-order forward-dif ference formula and 
Equation 6.2 is a second-order forward-dif ference formula. Both first- 
and second-order backward-difference formulas or alternating forward- 
backward formulas could also be tried. 

Second, more grid points could improve the resolution of the shock 
layer detail. By increasing the number of grid points in the radial 
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direction, better boundary layer resolution could be achieved making 
an accurate solution at higher Reynolds number possible. If the number 
of grid points in the cross-flow direction were increased, better reso- 
lution of the transmitted shock could be achieved. Because of the 
storage- saving procedure used in this study, the use of more grid points 
is possible and could be very beneficial. In addition, the resolution 
of the transmitted shock could be greatly improved by clustering the 
grid points in the cross-flow direction near the intersection line. 

Third, a better comparison with experiment could be obtained if 
the inflow boundary condition in the cross-flow direction could be made 
to more adequately match the experiment. This could be achieved by 
using a different inflow boundary condition or by comparing with an 
experiment in which the impinging shock strikes the bow shock suffi- 
ciently far downstream from the blunted end (about six diameters). 

Fourth, different shock interference patterns (Types IV and VI) 
ciuld be computed by varying the sweep angle (X). The Type IV inter- 
ference pattern might be extremely difficult to compute because of the 
subsonic cross -flow Mach number. 

Finally, other problems related to the present problem are given 
by the following: 1) Shock impingement with real gas effects, including 

both chemical equilibrium and chemical nonequilibrium. 2) Shock im- 
pingement where the body is a hemisphere. 3) A wing leading edge solu- 
tion in which the body is more accurately modeled as a wing (i.n., '.Tri- 
able body radius, taper, angle nf attack, n':c.). 
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IX. APPENDIX A: REAL GAS SOLUTION PROCEDURE 

Presented in this appendix is the solution procedure (including 
the shock jump conditions) when real gas effects are included. The 
fluid is assumed to be air in chemical equilibrium. The solution pro- 
cedure is identical to the procedure used for a perfect gas case v'th 
the following two exceptions: First, all thermodynamic properties of 

the fluid are determined from curve fits designed especially for use 
with finite-difference methods [74,751. The curve fits available are 
indicated by 

P = P(e, p) 
a - a(e, p) 

T = T(e, p) A. 1 

h = h(p, p) 

T = T(p, p) 

where p is the pressure, e is the internal energy, p is the density, 

a is the speed of sound, T is the temperature, and h is the static 

enthalpy. These curve fits are valid for temperatures up to 25,000 °K. 

-7 3 

and densities from 10 to 10 amagats. 

Second, the expressions for p^ and in the shock jump conditions 
(Equations 3.14 and 3.15) are not valid when the fluid is assumed to 
be in chemical equilibrium across the bow shock. Closed-form analyti- 
cal expressions for and for fiu: real gas case cannot be found 
because a closed-form analytical expression for the equation of state 
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does not exist. Therefore, for the real gas case, an iterative proce- 
dure is required. 

The freestream density (p w )j the freestream pressure (p^ , and 
the freestream static enthalpy (h ) are all known and held constant, 

DO 

just as in the perfect gas case. The pressure just downstream of the 

bow shock (p£) is computed by exactly the same technique used in the 

perfect gas case. The iteration procedure, which is used to determine 

p£ and V N , begins at this point. The density ratio (DR = p^/p^) is 

the basic iteration parameter. Both lower and upper bounds for the 

density ratio (DR . and DR ) must be computed at the start of the 
mm max 

iteration and are given by 


DR . = 0.0 

mm 


DR 


max 


h 2 Y - 1 
— “T+ 1 

P« Y +1 
2z± + Z 2 

Y +1 


A. 2 


A. 3 


Notice that the upper bound density ratio (Equation A. 3) is obtained 
from the perfect gas formula. 

Next, the exact value of the density ratio (DR) is approximated by 

DR = (DR . + DR )/ 2 A. 4 

mm max 


Using this density ratio, a value for is computed from 


P 2 " P * 


V N \l p (1 + DR) 


A. 5 
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which is an exact shock jump relation. Next, a value for is 
computed from the curve fit formula 


^2 ^ 2 ^ 


A. 6 


With these values of V N and h 2J the density ratio is recomputed from 


/^l _tl 2 

DR = 2 .. ■ + 1 

com \l „ 2 


A. 7 


V. 


N 


which is another independent relation obtained from the shock jump 
conditions. At this point a convergence test, which involves a com- 
parison between |DR - D ^ com S ant * a predetermined error tolerance, is 
made. If the test is passed, the iteration is stopped, and the current 
values of p„ and V„ are used to compute the velocity components (u ? , 
u _ , and u „) and the radial shock velocity (r s ) from Equations 

3.16-3.19. If the test fails, either DR . or DR is updated 

min max 

according to the following scheme: 


If 

DR 

com 

> DR 

then 

DR . - DR 

mm 

A. 8 

If 

DR 

com 

< DR 

then 

DR = DR 

max 

A. 9 


With the newly computed value for either DR . or DR the process 
J r mm max 

given by Equations A.4-A.9 is repeated until the convergence test is 
satisfied. 

During the iteration process a special problem associated with 
Equation A. 7 may occur. If the density ratio given by Equation A. 4 is 
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far enough from the desired result the radical in Equation A. 7 will 

be negative. This situation is alleviated by returning to Equation A. 4 

with DR set equal to the current value of DR. 
max n 

The process just presented is known as a bisection type of itera- 
tion. Convergence is guaranteed providing that the final value for DR 

lies between the initial values of DR . and DR . Because this rou- 

mm max 

tine only halves the solution interval for each iterations many itera- 
tions are required to obtain reasonable accuracy. To reduce the number 
of iterations a linear interpolating procedure for obtaining improved 
density ratio approximations has been added. A schematic of this pro- 
cedure is shown in Figure A1 where the density ratio values lie along 
the horizontal axis and the errors (ER = DR - DR ) lie along the 
vertical axis. Notice that both the minimum and maximum density 
ratios and the corresponding errors for these density ratios are indi- 
cated in Figure Al, Improved approximations for DR are obtained by 
replacing Equation A. 4 with 


DR = DR . 

mm 


4- ER 


DR - DR . 
max mm 


min ER . + ER 

mm max 


A. 10 


where ER * and ER are determined from 
mm max 


SR 


min 


= DR 


min 


- DR 


com 


A. 11 


ER 

max 


DR 


max 


-DR 

com 


A. 12 


The use of this linear interpolating scheme may speed up the iteration 
process by as much as three times. 




Figure Al. Lineai interpolating procedure used in the shock jump iteration 
scheme. 
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The zero normal pressure gradient boundary condition which was 
used to compute wall values of density if given by 

P NI = P NI-i e NI-l /e NI A ’ 13 

where the subscripts NI and NI *• 1 represent the values at the wall and 
one grid row above the wall, respectively. This boundary condition was 
used in the real gas computation, even though its formulation requires 
a i 'feet gas assumption. For physically realistic cold wall solutions 
in the real gas range, however, the temperature in the vicinity of the 
wall (i.e,, in the lower part of the boundary layer) is probably low 
enough to make this boundary condition valid. 

Sutherland's viscosity law, in conjunction with a constant Prandtl 
number assumption, were us; d to compute the coefficients of viscosity 
(p) and thermal conductivity (k) . These assumptions are generally not 
valid for real gas calculations. However, other alternatives, such as 
curve fit approximations (for p and k) , are not available at the present 
time . 

The above real gas procedure was used in conjunction with the 
previously presented finite-difference procedure to compute a sample 
infinite cylinder solution. Because of the inaccurate method for 
computing the real gas transport properties, the solution computed with 
the real gas procedure was chosen to be in the perfect gas range. Com- 
parisons made between the real gas solution and the corresponding per- 
fect gas solution were very good. The real gas computation required 
50 percent more execution time on a per time step per grid point basis 
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than the perfect gas computation. By this means the real gas procedure 
was tested. Before any effective real gas computations can be made 
with this procedure, a fast approximate method for evaluating the 
transport properties of air in the real gas regime must be found. 
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X. APPENDIX B: SHOCK JUMP CONDITIONS 

The derivations of the three-dimensional shock jump conditions as 
well as the three-dimensional radial .shock velocity equation are pre- 
sented in this appendix. The resulting equations are in a form con- 
venient for computer application. Figure B1 shows a sketch of the 

local shock-fixed coordinate system (1 ,1 ,1 ) along with the standard 

N J. D 

cylindrical coordinate system (i ,1 ,1 )• 

it y z 


A. Three-Dimensional Shock Jump Conditions 


The resultant velocity vector on the downstream side of the bow 
shock (V 2 ) may be expressed as the vector sum of 'he freestream veloc- 

1^ ■ i 

ity vector (V ) and some arbitrary velocity vector increment (AV) as 

CO 

given in Equation B.l. 



V + AV 

m, S S 


B.l 


Equation B.l is written in the shock-fixed coordinate system (see 

Figure Bl) . V and V. are given by 
co. s Z « S 


V 

00, S 


A . 
N , co N 


A , 

u x.^i + 


A 

U_ 

B, to B 


B. 2 



U N,2 1 N + U T, 2 1 T + u B,2 l B 


B. 3 


where u M , u_ , and u_ are the components of V in the shock 
N, 00 T, CO B, TO CO, s 

normal, tangent, and binomial directions, respectively, and u u , 

IN J Z 1 J L 
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Figure Bl. Shock-fixed and body-fixed coordinate sysLems. 
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and u are the components of V in the shock normal, tangent and bi 
B , 2 2 , s 

normal directions, respectively. 

Unon substitution of Equations B.2 and B.3 into Equation B.l and 
applying the fact that tangential components of velocity pass through 
shock waves unaltered yields 




U N,« 



B.4 


Note that the arbitrary velocity vector increment (AV g ) is normal to 
the bow shock as expected. 

The relative velocity equations for the and velocity vectors 
between the shock-fixed and the body-fixed coordinate systems can be 
written as 




B.5 


V 

CO, S 




B. 6 


where is the relative velocity of the bow shock with respect to 

the body and V , and V are given by 

2 « D 



A 

U 

r,2 r 


+ u 


A , 

, 2 1 0 + 


A 

U z,2 X a 


B. 7 


V 


>,b 


A A , A 

1 + U J- n +U 1 

ror 0, <» 9 2,0 z 


B. 8 


where u . , u„ n , and u „ are the components of V* , in the body coor- 

X" y 2 H y 2 Z j / 2 ^ D 

dinate system and u , u„ , and u are the components of V .in 

y r , 0, ca’ Z , » c CO, b 

the body coordinate system. Substitution of Equations B.4, B.S 
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and B.6 into Equation B.l yields 


A 


V 0 , = V . + (u T _ - u )i' 
2,b co,b N,2 N,»N 


B.9 


By taking components of Equation B.9 in the r- s 0- 3 and 2 -directions 
yields 


r,2 

= v • i = 
' 2,b r 

u + 

IT j 3 

CM 

s 

3 
*' — ' 

“ U N, « 

0,2 

Ti A 

2,b 

u n + 

0, “ 

(U N,2 

" U N,® 

2,2 

7? A 

V 2 ,b’ 1 z 

u + 

Z, CD 

(u H,2 

U N, GO 


,A A . 


, A A 


, A A 


B. 10 


B. 11 


B. 12 


The three components u , u and u „ are the desired velocity 

^ j Z 1 j Z Z j Z 

components just downstream of the bow shock in the body coordinate 
system. More useable expressions for the second terms on the right 
hand sides of these three equations must now be found. 

The general expression for the bow shock surface is given by 


r g = f ( 0) z) B.13 


or written in a different way 


F(r, q,z) = r g - f ( 0,z) =0 


An outward unit normal vector to the shock surface can be found by 
using 


A VF 

l.. - “ 


B. 14 


l N | VF | 


B. 15 
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whore the v() operator in cylindrical, coordinates is defined by 


,o = &li + i|ili a + dli 

w flr r r^f) ft 3 z z 


B. 16 


which yields 


A 

i 


( r sA)i D - r 


A 


s fl' s^Q * s z z 


N 


B. 17 


fl + (r S0 /r s ) + (r Sz )‘ 


Taking components of Equation B.17 in the r-, fl-, and z-directions 
yields 


A A 
1 * 1 
N r 


1 1 + (r s >J + 


e s' 



B. 18 


A A 
1*1 = 

N fl 


lx 

fl S 


^ + (r s fi /r s ) + < r s 2 ) 


A A 
1 * 1 = 
N z 


-r, 


1 + < r s 0 /r s > + <r s 2 >' 


B. 19 


B. 20 


Next, an expression for (u„ „ ’U„ ) in terms of freestream quan- 

N,2 N s o=> 

titles and the pressure just downstream of the bow shock (Pg) is de- 
sired. To this end the Rankin e-Bugon lot equations are used and are 
given by [58] 


P m + 
00 


P U VT 
r CO N, CD 


P 2 + P 


2 U N,2 


B. 21 



Combination of these two equations yields 


u„ „ - u 


P»-P 2 


N,2 N,oo p u 


B. 23 


N, co 


Substitution of Equations B.18, B.19, B.20, and B.23 into Equations 
B.10, B.ll, and B.12 yields the shock jump conditions for the velocity 
components in body-fixed coordinates, which are given by 


u _ = u + 


p .- p 2 


r ’ 2 r ’ c ’ P« U N. 


1 1 + < r « e /r s ) + (r s 2 >‘ 


B. 24 


u „ „ = u „ + 


p .- p 2 


■ r s> 


e s 


9>2 R,co p u , ? 9 

“ 1 1 + (r 8 Jr ) 2 + (r s ) 2 


Pi s' 


B. 25 


u + 


P P, 
cc Z 


-r< 


z , 2 a , ® p u 
* H CO Nj CO 


1 + (rs fl /r s )2 + (r s z )2 


B. 26 


To complete this set of shock jump conditions expressions for u 
and p£ in terms of freestream quantities and p^ are needed. The ex- 
pression for u^ m is usually written in terms of another quantity 
defined by 


N, CO 


V - v ■ (-i ) = 
n o,s k r 


“U 


N,® 


B.27 


For the case of a real gas where the flow is everywhere in chemical 
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equilibrium, including across the bow shock, closed-form expressions 
for and do not exist. Appendix A deals with this case by uti- 

lizing curve fits for the properties of air and an iteration procedure 
to compute and 

For the case of a perfect gas, closed-form expressions for V^, and 
are available [*54,64]. By using these expressions and the new 
definition for V , the entire set of three-dimensional shock jump 
conditions can be stated as 



B.32 
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B. Three-Dimensional Radial 
Shock Velocity Equation 


Figure B2 shows a sketch of a segment of the bow shock at time t 
and t + At and some of the key notation used in this derivation. The 
relative velocity equation for the freestream velocity vector (VJ 
between the shock-fixed and the body-fixed coordinate systems is given 
by Equation B.6. The freestream velocity component normal to the bow 
shock, measured in shock-fixed coordinates, and positive inward (V^) 
can be written with the aid of Equation B.6 as 


V = V 
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CO, s 


, * \ 
<-V * 


(V. 


>,b 


- w 


/ A \ 

( "V 


B. 33 


Equation B.33 may be simplified by using Equation B.8, B.18, B.19, and 
B.20 to yield 


-u 


+ u 


.(r Sn /r ) + u 


V N = 


r ’ co 8'” a 0 5' Z,cO 13 Z ^ 
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1 + <rs e /r s ) 2 ♦ (r..) 2 
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B. 34 


At this point reference to Figure B2 is necessary to obtain a more 
useful relationship for the last term of Equation B.34. Note that 
is the component of r s t ^ r t ^ ie shock normal direction which yields 

V /u - (r l • x ) x B. 35 

s/b s t r N N 


where r Sfc is the desired radial shock velocity. Substituting this 
result into Equation B.34 and simplifying with the aid of Equation B.18 


yields 



/. 



Figure B2. Notation used in the derivation of the three-dimensional 
radial shock velocity equation. 
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-u + u n (r g 7r ) + u r„ + r s 
_ r, oa R,co v s 6 s' z,» s z s t 

w I " 2 : 2 

fl+(r S0 /r s ) + (r Sz ) 


B. 36 


Solving this expression for r S{ _ yields the final result which is given 

by 


r s „ = I 1 + <*-7* ) 2 + < r s ) 2 + u 
s t N / s 0 s s z r,co 
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